% Sample H (highest) and SH (second-highest) values
% Unconditional on reserve price
% For a grid of n's
% Output matrices of [maxn, sims] Highest values and Second highest values
% Outputs values in HOMOGENIZED SPACE (g(z)=0)

function[HV,SHV]   = sampleHSHvalues(PDFb,CDFb,PAR)

     adjust             = 1/(integral(PDFb,PAR.minv,PAR.maxv));
     highfv             = @(s,n)     n.*(CDFb(s).^(n-1)).*PDFb(s);
     sechighgivenhighfv = @(t,s,n)   (n-1).*(CDFb(t).^(n-2)).*PDFb(t)./CDFb(s).^(n-1);
     highFv             = @(v1,n)    integral(@(s) highfv(s,n),PAR.minv,v1)*adjust;       
     sechighgivenhighFv = @(v2,v1,n) integral(@(t) sechighgivenhighfv(t,v1,n), PAR.minv,v2);  

     HV    = zeros(PAR.maxn,length(PAR.probsMC));
     SHV   = zeros(PAR.maxn,length(PAR.probsMC));
  
     for n = 1:PAR.maxn
         HVn = zeros(1,length(PAR.probsMC));
         for j = 1:length(PAR.probsMC)
            try
               HVn(j) = fminbnd(@(s) (highFv(s,n)-PAR.probsMC(j))^2,PAR.minv,PAR.maxv); 
            catch
               HVn(j)       = HVn(j-1);
            end
         end
         
         if n > 1
         SHVn = zeros(1,length(PAR.probsMC));
             for j = 1:length(PAR.probsMC)
                 try
                     SHVn(j) = fminbnd(@(t) (sechighgivenhighFv(t,HVn(j),n)-PAR.probsMC(j))^2,PAR.minv,HVn(j)); 
                 catch
                     SHVn(j)  = SHVn(j-1);
                 end
             end
         else
            SHVn(:)   = PAR.minv; 
         end
         HV(n,:)  = HVn;
         SHV(n,:) = SHVn;
     end
        % Impute any missing values - and transform to lnV
        HV  = exp(knnimpute(HV));
        SHV = exp(knnimpute(SHV)); 
end

 